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A  FINITE  DIFFERENCE  SCHEME  FOK  THE  NEUMANN  AND  THE  DIRICHLET  PROBLEM 

by 
K.  0.  Frledrlchs 

A  specific  feature  of  the  finite  difference  scheme  here  proposed 
Is  that  In  It  the  boundary  condition  and  the  differential  equation  are 
treated  simultaneously.   It  thus  differs  essentially  from  the  schemes 
for  the  Neumann  problem  referred  to  by  Forsythe  and  Wasow  (in  Section 
20.10  of  their  book)  and  the  scheme  of  H.  Keller  (to  appear),  in  which 
a  special  construction  is  employed  -to  set  up  a  finite  difference 
boundary  condition.   Our  scheme  is  close  to  those  in  which  the  bound- 
ary condition  is  derived  as  a  natural  one  from  a  variational  principle 
for  the  finite  difference  solution  and  which  tlierefore  has  a  symmetric 
coefficient  matrix.   Such  a  scheme,  indicated  for  Laplace's  equation 
by  Forsythe  and  Wasow  (Section  20. 5),  was  developed  for  systems  of 
equations  by  G.  White  (to  appear). 

Our  scheme  will  also  result  from  a  variational  principle,  but 
from  the  variational  principle  of  the  original  differential  equation 
problem,  simply  by  using  the  Ritz  method  employing  piecewlse  linear 
approximation  functions.   Thus  a  rather  uniform  treatment  of  boundary 
condition  and  differential  equation  will  result.   The  coefficients 
of  the  substitute  boundary  condition  will  be  given  as  the  areas  of 
sections  of  triangles  cut  out  by  the  boundary;  they  are  therefore 
not  sensitive  to  variations  of  the  direction  of  the  boundary.' 

An  advantage  of  our  approach  is  that  the  mean  convergence  of 
the  solutions  of  the  difference  equations  (including  first  differ- 
ence quotients)  to  that  of  the  differential  equation  is  implied  by 
the  general  theory  and  requires  no  special  proof.   Whether  or  not 
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the  scheme  Is  useful  In  actual  computation  Is  still  to  be  seen. 

At  first  we  shall  describe  the  scheme  for  the  Neumann  problem 
of  Laplace's  equation;  then  we  shall  discuss  its  extension  to 
systems  of  second  order.   In  two  appendices  we  shall  give  a  lower 
estimate  for  the  eigenvalues  of  the  matrix  involved  and  shall 
indicate  how  the  scheme  may  be  extended  to  the  Dirichlet  problem. 

1 .   The  Triangulation  Scheme  for  the  Neumann  Problem 

The  Neumann  differential  equation  problem  is  that  of  finding  a 
solution  <l)(x,s)  =  u(x,s)  of  the  equation 

Au  =  -f(x,s)  (1) 

;ion  i\_  which  satisfies  the  condition 


m  a  rei 


3u 


^--els)  (2) 

at  the  boundary  (S  ,    where  n  refers  to  the  interior  normal;  the 
given  functions  f(x,s)  and  g(s)  should  satisfy  the  condition 


//: 


f  dxdy  +  /   gds  =0,  (5) 

necessary  if  there  should  be  a  solution  of  the  problem. 
Under  appropriate  conditions  onWi_,  6  ,    f,    and  g  -  which  we  do  not 
formulate  here  -  a  solution  exists,  which  is  unique  -  if  a  condi- 
tion such  as 

//   udxdy  =0  (^) 

is    imposed.      This    solution  u  minimizes    the   expression 

I[<t>]    =   i  D[c^]    -    if     fO    dxdy    -    /      g<tds 

^^^  6 
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among  all  continuous  functions    ^    in    (/C  +   C>     which  possess  piecewise  con- 
tinuous first  derivatives.     Here    D[^]  =  D[^,^]    and 

D[u,  ^]  =    /    /     \  {  —  )     +(-r^)     (  dxdy.       Since  the  first  variation  vanishes  for 
J  J^    (     Sx  ay       ) 

all  functions    ^,     i.  e. 


the  relations 


and 


D[u.(<.]    =      /   /       f^dxdy    +      /       g^ds, 
D[u]    -     I  I        fudxdy    +      /       gu  ds 


I[u]    =      -  \  D[u] 


hold. 

Instead  of  the  Neumann  problem  proper  that  we  have  described  we  may 
just  as  well  consider  a  more  general  problem  involving  a  self-adjoint  differ- 
ential equation    Lu  =  f    with  a  boundary  condition    Mu  =  g    where    u,    f,     and    g 
stand  for  systems  of  functions  rather  than  for  single  functions.     The  operators 
L    and    M    should  be  associated  with  a  positive  definite  quadratic  form  in 
place  of  the  Dirichlet  integral    D[u].       Our  way  of  deriving  a  finite  difference 
problem  in  the  present  section  can  be  carried  over  to  this  general  case  and 
all  convergence  statements  to  be  made  are  valid  just  as  well.     A  specific 
discussion  of  the  finite  difference  problem  for  the  general  case  will  be  given 
in  Section  4. 

We  place  in  the    x,  y  plane  a  rectangular  net    N    in  the    x    and    y    di- 
rections   with  mesh  width   a    and    i3    and  pass  a  diagonal  in  each  rec- 
tangle.    Specifically,   we  pass  the  diagonals  according  to  the  alternating 
scheme  indicated  in  Figure  1.      To  each  net  point    P    we  assign  as 
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Its  star^p  the  rectangle  formed  by  the  triangles  of  which  It  Is  a 


vertex,  and  as  Its  reduced  star  -3  t,  the  intersection  of  the  star 


^  p  wl 


th  the  rei 


;lon  C^ 


t 
P 


have  points  In  common  w 


We  reject  all  vertices  whose  stars  do  not 
1th  the  (open)  region  (l{^  .      The  net  points 


P  \ 


<\ 


Figures  1,  2 
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retained  form  the  set    (P);      the  number  of  these  points  will  be  called    k    . 
To  each  point    P    in  this  set  we  assign  an  (unknown)  value    ^    .      The  finite 
difference  equation  to  be  set  up  will  be  an  equation  for  these  unknowns. 

To  set  up  such  an  equation  we  assign  to  the  set  of  values    (<S„)    the 

N 
function    ^    (x,  y)     which  is  linear  in  each  triangle  of  the  net  and  assumes  the 

values    ^       at  each  point    P    of  the  set    (P),     while  it  vanishes  at  every  re- 
jected vertex. 

We  denote  by    ^^(x,  y)    the  function  which  is  linear  in  each  such  triangle, 


assumes  the  value    1    at    (x,  y)  =  P    and  vanishes  at  the  boundary  of  the  star 

Q  N 

/O       and  outside  of  it.     The  function    ^    (x,  y)    assigned  to  the  set  of  values 


(i    )    is  then  simply  given  as 


N  r— 

^    (x,  y)     =     2_p  ^TD^P^^'^^' 


Evidently,   the  functions    rj    {x,y),     and  hence  the  piecewise  linear  function 

N  N  /O 

^    (x,  y),     are  continuous.     We  also  note  that    <^    (x,  y)  =  0    in    ^     implies 

^      =  0    for  each    P    in    (P).       We  now  consider  the  expression    I[§  ]    in  the  re- 

N  N 

stricted  class  of  functions    4       associated  with  our  net    N.       Evidently,     l[^     ] 

A 

may  be  regarded  as  a  quadratic  function  of  the    k       values    ^       generating  the 

N 
function    ^     .       This  functional  evidently  has  a  finite  lower  bound  since  obviously 

I[^^]    >     I[u]. 

N 
By    u       we  denote  a  function  for  which  this  minimum  is  attained.     It  follows 

that  this  functional    I((^)    has  a  minimum.     Evidently,    the  first  variation  of 

l[(^]    taken  for    (^  =  u       vanishes  for  all  functions    <f>     : 
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D[u^♦'']    =    IT    f/'  dxdy  +    f    sAs 


N  N 

For   4)      =   u      this   gives 

D[u^]    =    ff    fu^   dxdy  +    Tgu^ds 

and  hence 

iLu    J    =    -^-DLu    J  . 

The  solution  u   Is  unique  except  for  an  additive  constant.   For, 
clearly,  D[v  ]  =  0  for  the  difference  v   of  two  solutions  whence 
V  =  const.  In  all  triangles  which  Intersect   the  open  region  Kl  , 
and  every  admitted  point  P  Is  vertex  of  such  a  triangle.   The 
solution  u   can  therefore  be  made  unique  by  Imposing  a  condition 


such  as 


'<^ 


u  dxdy  =  0. 


The  strong  convergence  of  the  solutions  u   to  the  solution  u 
of  the  full  problem  Is  easily  established.   First  we  note  that  the 
function  u  can  be  approximated  by  functions  <t>      -   with  N  sufficiently 
large  such  that  D[4>  -u]  ,  //^  {<^   -u]   dxdy,  and  \^^\^   -u]   ds  are 
arbitrarily  small.   From  this  fact  it  follows  that  also  l[<l)  ]  can 
be  made  arbitrarily  close  to  l[u] .   This  is  then  true  even  more  so 
of  I[u^]  <  I[ct)^].   Thus 

I[u^]  — >l[u] 

if  the  net  N  is  sufficiently  refined.   Now  from  the  vanishing  of 
the  first  variations  for  u  we  find 
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D[u,u  ]    =  fu  dxdy  +  /   gu  ds  =  D[u  ], 


hence  D[u-u^]  =  D[u]  -  2D[u^]  +  D[u^] 

=  D[u]  -  D[u^] 
=  2l[u^]  -  2l[u]  -^  0 

if  the  net  N  Is  sufficiently  refined.   From  this  relation  one  derives 
the  strong  convergence  relation 


fj     (u-u^)^  dxdy  ->  0 


by  Polncare's  inequality  and 

f     (u-u^)^  ds  ^  0 

by  using  another  standard  inequality  (see  Courant^  Hllbert ,  Vol.  II, 
Ch.  VII) 0 

This  scheme  that  we  have  described  suffers  from  one  defect 
the  fact  that  there  are  two  types  of  net  points,  those  whose  star 


/O  -r,   consists  of  ei£ 


^^..^^^^^   w^  ..o-ght  and  those  whose  star  consists  of  four  tri- 
angles of  the  net.   Accordingly,  as  will  be  seen,  the  finite  differ- 
ence equation  which  will  result  from  our  scheme,  is  different  for 
different  types  of  points  P.   An  Improvement  could  be  obtained  In 
two  different  ways. 

One  may  compute  the  solutions  for  each  of  the  two  devices 
of  triangulation  and  then  take  the  mean  values.   The  difference  of 
the  two  solutions  may  be  used  to  estimate  the  accuracy  of  the 
solution  so  obtained. 
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Or  one  may  to  begin  with  employ  the  mean  value  of  the  two 
functions  Tip(x,y)  associated  with  the  two  ways  of  trlangulatlon.   In 
this  case  all  points  are  of  the  same  type,  but  the  finite  difference 
equations  are  more  complicated   (9  point  Instead  of  5  point 
equations) . 

After  having  made  either  of  these  Improvements  It  would  be 
possible  to  show  (e.g.  by  the  methods  of  Courant,  Frledrlchs,  Lewy) 
that  the  solutions  u   of  the  difference  equations  converge  In  any 
interior  region  polntwlse  to  the  solution  u  and  Its  derivatives  pro- 
vided the  function  f  satisfies  appropriate  differentiability 
conditions . 

We  shall  not  elaborate  on  these  matters  of  improvements 
(although  we  shall  refer  to  them  later  on) ;  rather  we  shall  confine 
ourselves  to  the  simple  scheme  as  described  and  discuss  in  detail 
the  nature  of  the  associated  finite  difference  equation. 

2.   The  Finite  Difference  Equation 

The  condition  that  the  first  variation  should  vanish  for  the 
solution  u  of  the  minimum  problem  is  equivalent  to  a  system  of  .; 
kj,  equations  for  the  k^^  values  Up.   This  system  may  be  Interpreted 
as  a  finite  difference  equation. 

To  describe  the  specific  nature  of  this  equation  we  set 

omitting  the  Index  N  here  and  in  the  following.   We  then  find 

D[<t)]  =  D[^p  0p  Tip] 

-  Up  dp  <fp+2];^Z:p.4p  dpp,  <J'p,  *p    > 
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where  we  have  used  the  abbreviations 

dp  =  D[Tip]    ,    dpp,  =  D[rip,,Tip]. 
Evidently,  dp,p  =  dpp,.   The  associated  bilinear  form  is 

D[$,(l)']  =   ^p  dp  (t-p  <t)p  +  2;>^J^p,^p  dpp,  cf'p,  0p. 
Hence  the  system  of  equations  can  be  written  as 

dp  Up  +  ^ p  , ^p   dpp I  Up ,  =  ep 


where 

er,  = 


^  =   [[      f  Hp  dxdy  +   j     g  Tip  ds. 


>''(R     ^         ''© 


The  coefficients  of  this  equation  are  easily  evaluated.   In 
doing  this  we  shall  find  that  the  coefficients  dpp,  vanish  unless 
the  point  P'  is  a  nearest  neighbor  of  the  point  P,  i.e.  if  it  differs 
from  this  point  either  only  by  ±  a  in  its  x-coordlnate  or  only  by 
±  P  in  its  y  coordinate.   Denoting  these  points  P'  by  P+^j  P+g  and 
the  corresponding  coefficients  dpp,,  by  d^  ,  d_^g  we  therefore  may 
write  our  equation  -  omitting  the  subscript  P  -  as 

du+du  +d  u   +  doUQ  +  d  o^  o  =  e 
a  a    -a  -a    p  p     -p  -P 

for  each  point  P  of  the  set  (P).   Moreover  we  shall  find  that 

In  evaluating  the  coefficients  we  must  distinguish  two  types  of 
points  P.   We  say  the  point  P  is  of  Type  I  if  its  star^p  consists 
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of  eight  triangles  as  indicated  in  the  figure 


Figure  3-i 
(Star/^p  of  Type  I) 

The  point  P  v;ill  be  called  of  Type  II  if  its  star  2>p  consists  of 


four  triangles  as  in  the  next  figure. 


Figure  Jg 
(Star^p  of  Type  II) 

The  reduced  stars  jSp  then  consist  of  eight  and  four  possibly 

mutilated  or  absent  triangles.   Possible  cases  are  shown  in  the 

figure  below. \ 


,      ~-      Figure  4 
Possible  star  of  /5  p  of  Type  I  Possible  star 


/p  of  Type  II 
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In  each  of  these  two  cases  we  denote  by  A_^   the  two  possibly 
mutilated  or  absent  double  triangles  with  altitude  In  the  x- 
dlrectlon  as  shown  In  the  figures  below: 

\ 


a  -a 

Figure  5 


Similarly,  we  denote  by  A  _  the  two  possibly  mutilated  or  absent 

double  triangles  with  altitude  In  the  y-dlrectlon: 

I 
\ 

V 


We  now  maintain  that  In  both  cases 


d 


-2 


-2 


±a 


■^    ^±a    '         ^±P  ==  -^    ^±P  ' 


and 


d  =  a  ^(A^  +  A_^)  +  p  ^(A^  +  A_p)  , 


so  that  our  equation  can  be  written  as 
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-2  -2 

a"   A  (u-u  )  +  a"   A   (u-u   ) 
a    a         -a    -a 

{*) 

+  ^"^  A^(u-u  )  +  p"^  A_g(U-U_g) 


To  prove  these  statements  we  observe  that  In  case  P  Is  of  Type  I 

the  function  t]  =  Tip(x,y)  Is  Independent  of  y  In  A^and  Independent  of 

X  In  A_^Q ,  while  b  t]   =   cl~      In  A_^   and  b   r\   =  ^~      In  A_^q  .   The  formula 
±p         X  '  ±a      y  ±P 

for  d  is  then  an  Immediate  consequence.   Next  we  note  that  the  prod- 
uct  d  T]„S  ri^,  differs  from  zero  only  If  P '  =  P_^   and  then 
X  'P  X  'P '  ^  ±a 

o 
^x^P'  "  ~^x^P  ^hil®  ^y^P'  "  \^P  "  °'   ^s^ce  d^^  =  a"  A^^  results. 

Similarly,  d^g  "^  ^'^   A^o- 

In  case  P  is  of  Type  II  we  find  (3^il  )^  =  a"^  and  (S  np)^  =  ^"^ 
in  each  of  the  four  "triangles"  forming  the  reduced  star  /j  p.   The 
formulas  for  d  then  follow  Immediately.   Furthermore  we  have 

_o 

^  Tc^  ^E.1  =  -ct   and  d  ri^S  Tr,  i  =  0  for  P '  =  P^   so  that 
X  'P  X  'P  '  y  'P  y  'P  '  ±a 

-2  -2 

d,   =  -a   A,  .   Similarly,  we  find  d^^  =  -P   ^^a-      Thus  it  is  seen 
±a         ±a  '^  ±^         ±p 

that  the  formulas  given  above  are  correct. 

From  the  derivation  of  our  finite  difference  equations  it  is 
clear  that  they  have  a  solution.   As  was  stated  earlier  the  solution 
is  unique  except  for  an  additive  constant;  clearly  u  =  constant  in 
a  solution  of  the  homogeneous  equations.   The  dependence  relation 
of  the  left  member  is  obtained  by  adding  all  of  them,  as  is  readily- 
seen.   The  right  members  satisfy  automatically  the  corresponding 


relation  2_p®p  ^   ^   since  this  relation  is  nothing  but 

//   fdxdy  +  /  gds  =  0  by  virtue  of  ripTip(x,y)  =  1  for  x,y  in  (^  . 
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To  solve  our  equations  numerically  one  may  omit  one  of  them 
and  choose  one  of  the  values  Up  arbitrarily.   Instead,  one  may 
perhaps  follow  the  "advice  to  the  programmer"  given  by  Forsythe 
and  Wasow  (Section  25-95  P-  576). 

It  is  evident  from  the  derivation  that  our  finite  difference 
equations  are  of  positive  type. 

5.   Interpretation  of  the  Finite  Difference  Equations 

To  Interpret  the  relationship  between  the  form  {*)    of  the 
finite  difference  equation  and  the  differential  equation  we 
consider  first  the  case  that  the  point  P  is  strictly  in  the 
Interior;  by  this  we  mean  that  the  whole  star  />  p   is  in  the  closed 
region  (\  +  Oj>  so  that  /ip  =  /Jp.   In  this  case  the  triangles  A_^ 
and  A_|_Q  evidently  have  the  area 

ha  =   ^±p  =  «^ 
so  that  equation  (*),  after  dividing  it  by  -a$ ,  takes  the  form 

a~^[u  -2u+u   ]  +  p"^[ua-2u+u  «]  =  -(ap)~  e. 
a      -a    *^     P      -p 

The  left  member  here  corresponds  to  the  Laplaclan  of  u.   However, 
the  right-hand  side 

(ap)"^e  =  (ctP)'^j(L  Tip  f  dxdy 
does  not  correspond  to  the  right  member  of  the  differential  equation. 
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Clearly,  the  ratio  // C?^p  f  dxdy  /  //^  rip  dxdy  Is  a  mean  value  of 

f  and  hence  corresponds  to  fpj  but,  since  the  area  of  a  star  of 
Types  I  and  II  is  4aP  and  2aP  respectively,  one  finds 


JJ     rip  dxdt  -  (4/3)ap 
'  /p 


in  case  P  is  of  Type  I, 


=  (2/3)ap        in  case  P  is  of  Type  II. 

The  right  member  (a^)"   ep  thus  corresponds  to  (4/3)f  and  (2/3)f 
respectively. 

This  fact  does  not  contradict  the  statement  made  above  that 
the  solutions  of  the  finite  difference  equations  together  with 
first  difference  quotients  converge  strongly  to  the  solution  and 
its  first  derivatives;  but  it  shows  that  the  second  difference 
quotients  in  general  do  not  converge  strongly  to  the  second  derivatives  of  the 
solution.  . 

It  is  easy,  though,  to  eliminate  this  shortcoming  by  modifying 
our  procedure.   In  fact,  either  of  the  two  improvements  mentioned 
at  the  end  of  Section  1  would  serve  this  purpose.   We  describe  in 
more  detail  the  first  one.   In  setting  up  our  procedure  we  have 
drawn  the  diagonal  of  one  of  the  rectangles  of  sides  a,  ^  in  an 
arbitrary  way.   We  could  have  drawn  it  the  other  way.   Points  of 
the  first  type  would  then  become  points  of  the  second  type  and 
vice  versa.   Let  u  and  u   be  the  two  solutions  of  the  two  finite 
difference  equations  so  obtained.   Since  the  left  members  of  the 
equations  satisfied  by  u  and  u   are  the  same  for  strictly  in- 
terior points  the  function 
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u 


(1/2)  u^  +  (1/2)  u-^-^ 


will  satisfy  at  such  points  an  equation  with  the  same  left-hand  side 
and  the  right-hand  side 

(l/2)e^  +(l/2)e^^, 

which  evidently  corresponds  to  aP  f  as  desired. 

At  present,  however,  we  continue  to  discuss  the  nature  of  the 
finite  difference  equation  resulting  from  one  of  the  two  procedures 
of  triangulation. 

Suppose  the  point  P  is  not  a  strictly  interior  point.   Then 
the  finite  difference  equation  does  not  correspond  to  the  differen- 
tial equation;  it  corresponds  to  a  linear  combination  of  differential 
equation  and  boundary  condition  if  the  point  P  lies  inside  the  region 
J^C^  and  just  to  the  boundary  condition,  except  for  a  factor,  if  it 
lies  outside  -^C   Here  we  assume  that  the  boundary  has  no  corners. 

In  fact,  in  discussing  these  matters  we  shall  even  assume  that 
the  boundary  is  straight  in  the  neighborhood  of  the  point  P.   Within 
terms  of  the  next  order  in  the  mesh  width  the  results  remain  true, 
if  the  boundary  is  curved  provided  its  normal  is  continuous  and  the 
mesh  widths  are  chosen  sufficiently  small.   The  two  components  of 
the  interior  unit  normal  will  be  denoted  by  i   and  m.   The  boundary 
conditions  then  can  be  written  as 

-2Su  +  mSu=-g. 

X      y 
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Within  terms  of  higher  order  we  may  replace  both  difference 
quotients  a"  (u^-u)  and  a"  (u-u  ^)  by  the  same  derivative  ^^u,  at 


a 


X 


a  representative  boundary  point,  and  p   (u^-u)  as  well  as  ^   (u-Ug) 
by  the  same  derivative  bu   there.   At  the  same  time  we  replace 
e  =  /   grids  +  /  /  _  f  rids  by  the  leading  term  /_  grids .   This  being 

done  the  finite  difference  equation  goes  over   (within  terms  of 
higher  order)   into  the  equation 

{**)  a'^ik    -A   )a  u  +  p"^(Ap.-A  p.)S,u  =  -/   grids, 

a   -ci  X         P   -p   y     "^/O 

Let  us  first  assume  that  the  point  P  is  in  the  exterior  of  the 
region  OC-   Let  us  further  denote  by  I  and  II   the  cases  in  which 
the  boundary  (D  intersects  two  adjacent  sides  of  the  star  ^  -p, 
and  by  I  and  II   the  cases  in  which  (q    intersects  opposite  sides 
of  ^ 


(I^  and  11^) 


Figure  7 


(A  -A  „) 
a        -/3 

(A„-A   ) 
p    -a 


(I^  and  11^) 
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We  now  maintain  that  In  each  of  these  four  cases  the  relation 


so  that  equation  {**)    can  be  written  In  the  form 


s 


where 


^  ^  u  +  m  S  u 
X       y 


g  =  /^  gilds  /r  nds 


B       /  ^(3 

Is  a  mean  value  of  g  associated  with  the  star/ip. 

The  verification  of  the  statement  (i:^)  Is  a  matter  of  elemen- 
tary geometry;  for,  the  regions  A_^  ,  A^q  are  triangles  or  trapezoids 
since  \C5  was  assumed  straight' In  /Sp.   Avoiding  tedious  case  distlnc-^ 
tlon  we  can  prove  the  statement  simply  by  Green's  formula 


/   (dTido  +  dTiS<J))  dxdy  =-/   rilid   +md)  <t>ds 
--(^   X   X     y   y  Jg     X      y 

where  11  =  11   and  «t>  is  any  potential  function^.   Choosing  <t»  -    iy   -   mx 

2 

and  (t  =  ^x  +  my  we  find 

m  a"^(A  -A   )  -  ^  p"^(A„-A  «)  =  0, 
a   -a  p   -p 


i  a~^(A^-A_^)  +  m  p"^  ( Ap^-A_f^ )  = 


i  a~^(A  -A   )  +  m  p  ^(A^-A  o  )  =  /   ^1  ds , 
a   -a  p   -p     ' 


and  thus  formula  (:p^)  . 
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As  a  result  we  have  shown  that  the  relation 

a"^[A  (u  -u)  +  A   (u-u   ) ] 
a  a       -a     -a 

+  p~^[Ag(ug-u)  +  A_g(u-u_p)]  +  e  =  0 

for  an  exterior  point  P  agrees  in  lowest  order  with  the  boundary 
condition 


[i^u  +  mSu  +  g]  /  ^  rids  =  0, 


6 

Formulas  (^)  from  which  this  fact  was  derived  are  also  valid 
for  interior  points.   They  become  empty  for  a  strictly  interior 
position  and  then  /   rids  =  0.   In  case  of  an  interior  point  P  near 
the  boundary  we  like  to  interpret  equations  (<fl* )  not  just  as  a 
substitute  for  the  boundary  condition  as  we  have  done  it  for  ex- 
terior points. 

Suppose  then  that  the  point  P  is  in  the  interior  of  the  region 


uv  but  not  strictly  so,  i.e.,  suppose  that  the  star  ^  p  does  not  lie 
completely  in  ^\  .    Then  let  us  denote  by  A^  ,  A_^g  the  areas  of  the 
triangles  associated  with  the  full  star  /O -p   and  set  A'   =  A    -  A_^  , 
A^Q=  A^   -  A_^  .   In  a  similar  way  we  set  e*  =  L    f  r\   dxdy  and 


e'  1=  e*  -  e  =  -/   g  ri  ds  +  /  /  ^   .\f  ri  dxdy.   The  finite  difference 
equation  can  then  be  written  in  the  form 

-  5a"^[A'  (u  -u)  +  A'  (u-u  )] 
<^      a  a       -a    -a      ~v. 

+  p~^[A'(ug-u)  +  A'  (u-u_p)]  +  e'^ 

+  a"^[A*(u  -u)  +  A^  (u-u  )] 
a  a       -a    -a  ^ 

+   ^^[A*(ug-u)  +  A*  (u-u_p)  +  e*  =  0. 
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The  first  expression  here  is  evidently  of  the  same  form  as  the  finite 
difference  equation  for  an  external  point.   The  fact  that  exterior 
and  interior  change  roles  is  balanced  by  the  minus  sign.   Hence  this 
term  is  a  substitute  for  the  boundary  condition.   The  second  expres- 
sion is  the  same  that  would  result  if  the  star  ^  p  were  completely 
in  the  interior;  hence  it  is  a  substitute  for  the  differential 
equation.   Thus  we  may  say  that  for  a  not  strictly  interior  point 
the  finite  difference  equation  corresponds  to  a  linear  combination 
of  boundary  condition  and  differential  equation.   The  right  member 
of  the  finite  difference  equation  then  also  becomes  of  the  order  of 
magnitude  1  in  the  mesh  width;  for,  the  term  ///P  f  'H  dxdy  being  an 
integral  over  the  reduced  star/^  p  Is  of  at  least  the  order  of 
magnitude  of  the  areas  A  +  A_  +  Ag  +  A_   covering  Jj  ^.   The  term 
Lq  g  ^1  ds  is  of  the  same  order  since  t]Q'    is  of  this  order  in  /^ -p, 
when  U -> '  is  the  length  of  the  intersection  of  \1>   with  /6p.   Thus, 
after  the  division  described,  the  diagonal  terms  of  the  matrix  of 
our  finite  difference  equation  equals  1,  the  off -diagonal  terms  are 
less  than  1,  and  the  right  member  is  of  the  order  1. 

4.   Difference  Equations  for  Systems  with  Non-Constant  Coefficients 

As  was  indicated  in  Section  1,  our  general  scheme  of  deriving 
finite  difference  equations  and  our  general  convergence  statements 
concerning  their  solutions  carry  over  literally  to  the  general  case 
of  self-adjoint  systems  of  equations  of  second  order  with  a  natural 
boundary  condition.   The  specific  nature  of  the  finite  difference 
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equations  resulting  In  this  general  case  requires  a  special  discus- 
sion. 

By  u,  f,  g  we  denote  systems  of  k  functions,  and  by  a,  b,  c 
we  denote  k  by  k  matrices  which  may  act  on  such  functions  and  which 
are  also  functions  of  x,    y.   We  assume  a  and  b  to  be  symmetric  and 
positive  definite,  while  c  Is  required  to  be  so  small  that  the 
Integrand  of  the  "Dlrlchlet  Integral" 

D[u]  =  // ^  [S  u-a  b   u  +   b   u-c  ^  u 
J  J  /?   X     X     X     y 

+Bu.c'Su+Su-bdu]  dxdy 
y     X     y     y 

Is  positive  definite  at  each  point  -  with  reference  to  ]^^u,  ^  u  ( 
as  2k  Independent  variables.   The  matrix  c'  here  Is  the  transpose 
of  c;  the  dot  Indicates  the  ordinary  Inner  product. 
The  variational  Integral  Is  then 

(l/2)D[u]  -  //  f-u  dxdy  -  /  g-uds. 

The  finite  difference  equation  Is  obtained  by  Inserting  the 
function 

and  varying  with  respect  to  Up.   The  result  Is  the  equation 


LpU  -  2 pi   pp I  ^p '  ~  ^p' 


where 


'  i/f^ '  ^p '^''^^  *  4  "^  "^p '''■ 
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while  the  matrices    K       ,    are  given  by 


Kpp,     =    Kpp,     +    Kpp,     +    K^p,     +    Kpp, 


with 


f-    r 


a 


Kpp.     =    -    /   /        a    r,p8^rip,c'dxdy, 

K^PP,    =   -yV     a^npB^np.bdxdy. 

Clearly,     K  vanishes  unless    P'    belongs  to  the  star  of    P.      In  general, 

therefore,   the  equation    L    u  =  e       is  a  nine- point  finite  difference  equation. 

In  order  to  describe  the  form  of  the  expression    L    u    more  specifi- 
cally we  introduce  the  points     P        =(+Q',0),     P    o  =  (0>  +/3),     P  ^=(+0-,  +3), 
where    P  =  (0,  0)    is  assumed.     Further  we  set 

Kpp  =  K,  Kpp        =K^^,  Kpp       =K^  K  =  K^^  ^ 

+  0-        -  1^        ~  +ar+3        -    - 

From  the  relations 

X  tt         -a  X     +p         ff+/3         -a +3 

which  are  immediately  verified,    we  deduce  the  relations 


K^    +    K^    +    K^       =    K^„    +    K^     „     +    K^        „ 
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Similarly  we  have 


and 


K^    +    K^    +    K^Q     =    K^        +    K^        o    +    K^       0     =    0 


c  c  c  c  c  c 

K+Kq+K„=K         +K         n    +   k        0=0, 
p  -/3  +a  +tt +/3  +a-/3 


a  -Of  +|3  a+p  -a+p 


Using  these  relations  we  may  write  the  differential  equation    L„u    in  the 
form 


o4-p'  ?^4-p'  P+P^ 

(#)  K  (u-u)+K      ^(u       -u)     +    K^'J  (u      „-u^) 

a  a  -a         -a  a+p      a+p        p 

T^a+c'  ,  ,,a+c'  ,  ,,a+c'  , 

-a+p     -a+p        p  a-p      a-p        -p  -a-p     -a-iS       -p 

+    Ko      (u„  -  u)    +    K   p   (u   Q  -  u)    +    K      p(u    ,o-u   / 
p  p  -p        -p  ff+p     Qf+p        a 

T^t)+c,  .  ,,b+c    ,  ,  ,,b+c    ,  , 

+    K      r,(u      o-u   )     +    K       ,o(u       ,Q-u      )     +    K         „(u        „- u      ), 
a-p     a-p        a  -a+p     -a+p         -a  -a-p     -a-p        -a 

where  we  have  set 


a  c'  a+c'  b  c  b+c 

K+K=K  and  K+K=K. 


In  order  to  see  in  which  way  this  expression  is  a  substitute  for  the  differ- 
ential expression    Lu    for  strictly  interior  points     P    we  use  the  quantities 

a=l  c-1 

K         ,  .  .  .  ,  K  obtained  by  setting    a  =  1^  ....  c=l    and  then  introduce  the  mean 

values 
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c  c 

a  ,,b  K        o  K        o 

-    _        K  -    _        K  -  ±a±^  -,  ±Qf±g 

^    "        a  =  l    '        ^    "        b  =  l    '         '^4<*+^     "     ^c  =  l     '  ''-Hff+0    "     ^c'  =  l    ' 

which  can  be  done  since  the  denominators  here  are  positive.     The  mean  values 

a      ,  a        Q  ,  b   o  ,  •  •  •     will  be  assigned  to  the  points    (  +  —a,  0),  (  +  —a,  +  —3),  .... 
+a       +01 +p       +p  —  2  —  2       —  2 

~  ~  ~c' 

In  addition  we  introduce  the  quantities    K        ^    by  restricting  the  integration  of 

c'  "^C 

K        over     +  y  >   y      and  the  quantities    K        q    t)y  restricting  the  integration  of 

c 
K   o    to  the  part  of  the  star  in  which     +  x  >  x   ,     where    (x   ,  y  )    are  the  coordi- 
+3  ^  -      -     o  o   ■^o 

nates  of  the  center    P:      then  we  introduce  the  mean  values    c'       o  =  K        ^/K        -, 

c  c=l  1  1 

and    c        o=K        ^/K        „,     assigned  to  the  point    (  +  — a,  +— jS).       Finally  we 

observe  the  relations 


-to  -w+i3  -har-/3  +/3  ff+jS  -Qf-i-/3 


~c'  =  l      ^    ~c'  =  l  -c  =  l     _    -c  =  l 


c'=l  c'=l  c=l  c=l 

-kjr+/3  -for-/3 '  ar+/3  -Qf+3  ' 

which  hold  for    T)  =  r]  ,     t]  =  r\    ,     and    r)  =  n  =  (l/2)r7    +  (1/2)17    ,     as  is  readily 
verified. 

Altogether  we  then  can  write  our  expression  in  the  form 
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(##)  <^^|a^(u^-u)-a_^(u-u_^)| 


+    K^   o  1  a      „(u    ,  o-  Un)  -  a      _  q(Uq-  u        p)  ? 

a=l  (  —  — 

+    K      Q  )  a      q(u      o- u    n)  -  a         p(u    p-U         r,) 
a-p  (     ar-p     a-p        -p  -a-fi     -p       -a-p 


+    K^"^  j  b^(u^- u)  -b_^(u-u_^)| 

+    K    ~o  <  b      o(u      „-  u   )  -  b      n(u   -  U      q)  > 
a+j3  (     a+/3    Qf+/3       or  a-i3     a       or-p    ) 

+    K    ~   o  <  b        o(u        Q-  u      )  -  b        p(u      -  u        n)  ( 
-ar+/3  (     -a+/3     -tt+i3       -or  -Qr-/3     -a       -or-p    ) 

+    k""'"^  \  (c'    p-c'     p)(u   -u)  j    +    k""'"!,  S  (c'    ^o-c'       J(u      -u) 
a+/3    {       a+jS        ff-p      a  )  -a+p  (        -a+p        -or-p       -or 

+    K^  ~     j  c'     o(u      o~  ^o)  ~  C     q(u      q- u    p)  f 
a+3    (     Qf+i3     a+jS       /3  a-p    a-p       -p    ) 

+     k"^  ~o  <   C'        n(u  p-  Up)   -    c'        p(u  Q-  u    p)  > 

-a+/3  (     -Qr+j3     -a+p       p  -a-p     -a-p       -p    ) 

+    K^~o  )  (c      p- c        p)(Up-u)|    +    k'^"p  >  (c      p- c        p)(up-u)> 
Qr+/3  (       Qr+/3        -or+jS       jS  )  Qr-/3  (       cir-i3        -a-p       -p  ) 

c  =  l  (  -  -  ) 

+    K      o  S  c      o(u      „-  u  )  -  c        p(u        p-  u     )  ^ 

a+/3  (     a+^     a+p       a  -a+P     -a+P       -a    ) 

c  =  l  (  -  -  ) 

+    K      pic      p(u      p- u  )  -  c        p(u        p- u      )  i 

a-p  (     a-p     a-p       a  -a-p     -a-p       -a    ) 


In  case  we  take    rj  =  H      we  find 


K^'^p    =    K^'"p    =( +  !)(  +  !)  if    P    is  of  Type  I, 

+a+p  +a+p  —       — 

=    0  if    P    is  of  Type  II. 
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Consequently,  as  the  mesh  width  shrinks  to  zero  the  expression 
(ap)    Lp  goes  In  the  limit  over  Into 

S      ab      u+S      bS      u+^25      c'S      u+2S      c'd      u.   ^ 

XX  yy     c.yx      xyj 

according  as  P  Is  of  Type  I  or  II.   While  thus  the  terms  In- 

2       2 
volvlng  B   and  S   are  represented  as  they  should  be,  this  is  not  so 

for  the  terms  Involving  S  S  . 

X   y 


Ind 


In  case  we  take  t]  =  'li  =  (1/2)11  +  (1/2)11   we  f 
^T^  =  0/4)a-^p,  K^=^  =  (1/8)  a-lp, 

<ll^  =   <a:l=    (1/2)(±1)(±1). 

If  the  mesh  width  .shrinks  to  zero  we  now  obtain  from  (ap)~   L  In 
the  limit  the  desired  expression 

S      ah     u+db^u+S   c'b   u  +   S    cS  y, 
XX  yy  yx  X      y"  ' 

It  Is  thus  Indicated  that  one  should  work  with  rj  rather  than 
T]  .   After  all,  once  c  7^  0  the  finite  dlffei'ence  equation  will 
Involve  nine  points  anyway;  no  substantial  saving  would  be  gained 
by  working  with  t]  . 

So  far  we  have  assumed  that  the  point  P  Is  strictly  Interior. 
Next  let  us  suppose  that  It  Is  exterior.   The  problem  then  arises 
In  which  way  the  finite  difference  equation  simulates  the  boundary 
condition 
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u 

X 


(ia  +  mc')^  u  +  (ic  +  mb)S  u  =  g. 

X  y 

As  in  Section  3  v;e  lump  together  all  difference  quotients 
v;hlch  correspond  to  the  same  differential  quotient,  obtaining 
the  expression 

-s      -s.       C  fira+c  '  ,  T^a+c  '  ,  ^a+c  '  \    (v3-+c  '  ,  r^-a+c  '  ,  x<ra+c  '  \  (  :n 
pa^u  +  qayU,=  a   ^(K^   +  K^^^  +  K^_^    )  -  (K_^^^+  K_^  +  ^-a-?^ C  ^ 

We  then  should  determine  the  ratio  of  the  coefficients  of  S  u  and 
S  u. 

y 

since  the  right  member  Is  e  =  e  =  (l/2)e   +  (l/2)e    It  Is  a 
substitute  of  (aP)~  f  according  to  the  results  of  Section  3-   We 
assume  a,b,c,c'  constant  and  restrict  ourselves  to  employing  t]  =  r). 

Now,  one  easily  verifies  the  rela,tlons 
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so  that 

p  =  (l/2)a-l(A^  -  A^^  .  Af  -  A^^)a  +  [1/2)^-^^1  -   A^^  +  A^^ 

q  =  (l/2)a-l(A^  -  A^^  .  A^  -  A^^)  c  +  [1/2)^-^^1  -  A^^ 

+  A^"^  -  A^^)b; 


thus  we  obtain  the  desired  formulas  p  =  0(ia  +  mc'),  q  =  0(ic  +  mb) 
according  to  the  results  of  Section  3  with 'e  =(l/2)  0"^  +(1/2)0"'""'", 
and  0=1    nds . 

Since  the  right  member  e  in  lowest  order  is  of  the  form 

/    r   I  r   II  ^ 

(1/2)  /  T]  gds  +(1/2)  /  T]     gds  it  agrees  in  lowest  order  with  0g, 
according  to  the  results  of  Section  3-   Thus  it  is  seen  that  the 
finite  difference  equation  for  an  exterior  point  in  lowest  order 
agrees  with  the  boundary  condition. 

Clearly  then,  as  in  Section' 3?  the  finite  difference  equations 
for  a  not  strictly  interior  point  may  approximately  be  regarded  as 
a  linear  combination  of  boundary  condition  and  differential  equation. 

Finally  v;e  make  a  few  remarks  concerning  the  question  whether 
or  not  our  finite  difference  expression  is  of  positive  type.   By 
this  we  mean,  in  extension  of  the  notion  introduced  by  Motkln  and 
Wasow,  that  each  of  the  coefficients  of  the  values  at  the  eight 
boundary  points  of  the  molecule  is  a  positive  definite  matrix.   This 
question  may  be  of  some  interest  although  a  maximum  principle  is 
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not  necessarily  Implied  by  positive  type. 

We  confine  ourselves  to  the  case  of  constant  a,b,c  and  consider 

only  the  case  of  a  strictly  Interior  point  P.   Then  the  coefficients 

—  2      -2 
of  u  and  u^  except  for  the  factor  (l/''4)aP,  are  3ct  a  -  p~  b, 

-2     -2 
3P   b  -  a  a.   The  condition  that  both  of  them  are  positive  de- 
finite may  possibly  be  achieved  by  choice  of  a  or  |i.   The  coeffi- 
cient of  u^  _,□,  except  for  the  factor  {l/>\)a^   Is 
±a±p 

(l/2)a~^a  +  (l/2)3"^b  ±  {±  )  (a(3)  "^  (  c+c  '  )  • 

This  condition  that  this  matrix  be  positive  definite  Is  satisfied 
If  |c  +  c'l  Is  sufficiently  small. 

Locally,  these  conditions  can  always  be  satisfied  by  proper 
choice  of  a  coordinate  transformation,  assuming,  of  course,  the 
equation  to  be  elliptic.   Whether  or  not   this  can  be  achieved 
globally  Is  another  question. 
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Appendix  I 

Remarks  about  Error  Estimates 

Concerning  the  accuracy  of  the  solution  of  the  finite  difference 
equation  one  should  first  observe  that  the  coefficients  of  some  of 
these  equations  may  be  rather  small,  namely  If  the  center  P  lies 
so  far  outside  the  boundary  (B  that  the  area  of  the  reduced  star 
^p  Is  very  small.   To  remove  this  smallness  one  need  only  divide  the 
equation  by  the  coefficient  of  u,  which  is 

d  =  a"2(A^,  +  A   )  +  3"^(Ap.  +  A  p.). 
a    -a         p    -p 

The  resulting  equation  will  be  called  the  "reduced"  finite  difference 

equation . 

It  is  clear  that  the  coefficients  of  the  unknown  u,  ,  u^^  in  this 

±ct   ±p 

reduced  equation  become  less  than  1.   The  right  member  of  the  reduced 
equation  -  we  maintain  -  will  be  of  the  order  of  magnitude  of  the 
square  of  the  mesh  width  for  strictly  Interior  points  and  of  the 
order  of  the  mesh  width  for  points  near  the  boundary.   Here  and  in 
the  following  we  assume  that  the  ratio  a/3  remains  fixed  if  the  "mesh 
width",  h  =  max(a,3)  tends  to  zero.   This  right  member  is  given  by 

e  =  /  /    T]  f  dxdy  +  /     "Hgds  . 
The  first  term  on  the  right  here  is  of  the  order 


ff      Tidxdy  ±  ff      dxdy  ^  h^d; 

J  J   g  I  J  J   O    I 
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the  second  member  is  of  the  order  /    rids  _^  hd,  as  seen  from  formula 
t>^'     on  page  20.   Since  this  second  member  is  absent  for  strictly- 
interior  points  our  contention  follows. 

The  question  naturally  arises  how  the  lowest  eigenvalue  of  the 
matrix  of  the  reduced  finite  difference  equation  behaves  as  the  mesh 
width  tends  to  zero.   We  recall  that  the  lowest  eigenvalue  of  the 
reduced  matrix  of  the  standard  finite  difference  equations  for  the 
Dirichlet  problem  tends  to  zero  as  the  square  of  the  mesh  width.   We 
maintain  that  the  same  is  true  of  the  lowest  eigenvalue  of  the 
present  reduced  matrix. 

Actually,  the  question  concerns  the  second  eigenvalue  of  this 
matrix  since  the  lowest  one  is  zero.   The  second  eigenvalue  A   can  be 
characterized  as  the  minimum  of  the  ratio 

D  [J^p  Up  \yy~^   dp  Up, 

for  all  functions  }  ^p  Up  t]  satisfying  the  condition  )  p  dp  Up  =  0, 
which  expresses  orthogonality  to  the  lowest  elgenfunction  Up  =  con- 
stant.  Our  statement  then  is  that  'h-^/a^   remains  bounded  as  a^  ^  0. 

Note  that  the  coefficient  dp  of  Up  in  the  unreduced  matrix  may 
be  very  small  if  the  point  P  lies  outside  of  the  boundary  IB;    but 
that  need  not  contradict  our  statement  since  the  coefficient  of 
Up  in  D['>~p  Up  rip]  is  Just  as  small. 

The  statement  that  A-,/aP  is  bounded  away  from  zero  independently 
of  the  mesh  width  is  closely  analogous  to  Poincare's   inequality  for 
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continuously  differentiable  functions,    which  we  formulate  in  a  somewhat 
extended  form 

/   /      pu^dxdy     <     C     /  /      [q  O   u)     +  q   (9    u)    ]  dx  dy 
^/   ^H  R        ""     ""  ^     ^ 

for  all    u    with      /   /     pudxdy  =  0,     and  an  appropriate  constant    C  >   0. 

R 

Here    p      q   ,     and    q      are  positive  functions  of    x    and    y    which  satisfy 
X  y 

certain  conditions  which  we  shall  describe.     We  assume  -  for  simplicity  - 
that  the  region    R    contains  in  its  interior  a  (closed)  connected  region    R' 
where    p    and    q    are  constant  and  that    R    can  be  covered  by  four  regions 
Q      ,     Q         possessing  translates    q'     ,     q'        in  the     +  x    and     +  y    direc- 

^+x      +y  ^  1'^      ly  ~  ~ 

tions  which  lie  in    R'.       Furthermore  we  assume  that    p    does  not  decrease 
along  the  segments  traversed  by  these  translations.     Also  we  assume  that 

p  <    cq   ,     <    cq      with  an  appropriate  constant     +  c.       The  regions  swept  out  by 

-        X      -        y  - 

the  translations  will  be  denoted  by    R      ,     R 

+  x         +y 


In  order  to  prove  the  extended  Poincare  inequality  we  first  rewrite 

] 
R 


,    ,  2 

/   /     pu   dxdy    in  the  form 


//      //      p(x     y  )p(x2,y2)[u(x^,y^)-u(x2,y2)]  dXj  dy^dx^dy^ 
R        R  

2     /   /      p(x,  y)dxdy 

I 
using  the  relation      /  /     pudxdy  =  0.       Then  we  proceed  essentially  as  in 

the  standard  proof  of  Poincare's  inequality  (see  Courant,   Hilbert,   Vol.   II, 
Ch.    VII). 

Let    (x'y'i    and    (x' ,  y' )    be  the  translates  of    (x   ,  y  )    and    (x   ,  y  ), 
then  we  first  estimate 
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To  estimate  the  first  contributions  let  us  assume  that  (x.  ,y-i  )  was 
translated  in  the  x-direction  so  that  y '  -   y-,  and  x-,  ^   x'  =  x.,  +  r. 
Then  we  have 

|u-^-uj2  <  [J  1  a^u(x,y^)dx]^  <  rj   ^  |  a^u(x^  ,y^)  |  ^dx. 
^1  ^1 

Integration  over  x-,  ,y-,  and  use  of  p(x-,  ,y.  )  ^p{x,y-,)  yields 


dx^dy^  2  W/    /    P^x^y^) l^x^^^'^1^  '   dxdx^dy^ 


^x  ^x  ^1 


2  r/ 


_  -Q  ,,   p(x,y^) |a^u(x,y^) I   dxdy^  ^  r^  c  ^x^^x^^   dxdy, 

X 

where  r   is  the  diameter  of  CR.   The  contribution  from  the  second 
o 

term  can  be  handled  in  a  similar  way;  in  fact  this  contribution 
is  given  by  the  original  Poincare   inequality  since  p  is  constant 
in  R '  .   Adding  all  contributions  we  find 

//  p^p^lu-^-u^l  dx-^dy^dx^dy^ 

<  l8r^  c  \'{    pdxdy  \\   [q^(a^u)^  +  q  (S  u)^]dxdy 

2 
whence  the  desired  inequality  with  C  =  9r  c. 

It  is  clear  that  a  similar  inequality  can  be  established  for 

the  quadratic  forms  associated  with  finite  difference  equations, 

provided  the  coefficients  p,q  ,q  satisfy  conditions  similar  to  i 

X  y 

those  required  of  the  functions  p,q^,q„- 

X   y 
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First  note  that  the  expression  D[J^pUprip]  can  be  written  In  the 
form 

X  y 

Here  M   stands  for  the  midpoints  between  two  x-adjacent  net  points 

P^   (with  x-coordlnates  differing  by  a  and  same  y-coordlnate) ;  V  u 

Is  the  difference  quotient  a"  [ u^  -Up   ],  and  q  =  (ap)~  A_   with 

ir    X  X  -i-Ct 

X    -X  + 

reference  to  P   as  center.   Clearly,  q  =  1  if  one  of  the  points 

IE  yC  Jv 

P^   is  strictly  interior.   The  terms  Involving  y  have  a  similar 
significance . 

Setting  Pp  =  (ap)~  dp  we  have 

dp  Up  =  ap  y"p  Pp  u. 


p  ^p  ^-p        ^y  /_ ^p  Pp  >^p  • 

-2     -2 
Clearly,  Pp   ^  2(a   +  ^   )  q^^,  for  P_^   strictly  in  the  interior. 

±x 
The  net  points  near  the  boundary  must  be  assigned  to  regions 

Q   and  Q   .   We  stipulate  that  an  interior  near  boundary  point  P 
±x      ±y 

is  admitted  to  Q  only  if  ^  >  0  in  its  star;  for  an  exterior- 
point  in  Q  in  addition  the  inequality  a   H   ^  ■;=p    |m|  should  hold. 
A  little  consideration  then  will  show  in  all  these  cases 

A  +  A    <  4A    ,   Ac  +  A  a  ""  ^A  • 
a    -a  —   a,      p    -p  —   a 

Denoting  by  M  the  midpoint  between  P  and  P  we  have  qjy,  =  (aP)~   A 
and  hence 

Pp  =  (ap)"^  dp  1  4{a~2  +  p~^)   qj^  =  cqj^. 


36 


This  Inequality  thus  holds  for  all  net  points. 

It  is  also  seen  that  the  value  of  p  does  not  Increase  as  a 

point  P  In  Q  moves  in  the  x-direction  Into  Its  translate.   Therefore 

/ 

our  derivation  of  the  extended  Polncare  inequality  can  be  carried 

over  to  the  present  case  for  the  domain  formed  by  all  (interior  and 
exterior  admitted)  net  points.   Let  r,  be  a  bound  for  the  diameter 

of  this  domain.   Then  the  resulting  Inequality  is 

2    -1      -1 
aP/A-,  ^   36r-,  (a~  ^  +  P~  a)  .   In  fact,  the  right  member  is  bounded 

in  the  mesh  width. 

The  method  of  estimating  the  eigenvalue  A-,  which  we  have 
described  can  immediately  be  carried  over  to  the  case  of  systems 
of  equations  provided  the  integrand  of  the  quadratic  form  D[u]  is 
positive  definite. 

In  the  case  of  the  equation  of  elasticity  without  boundary 
constraints  this  integrand  is  only  semidef inite ,  but  vanishes  only 
for  infinitesimal  rotations:  i.e.  for  u  =  03  X  x  where  cu  is  any 
constant  vector.   Instead  of  Polncare 's   inequality  one  must  use 
Korn's  inequality  (see  Payne  and  Weinberger  and  references  given 
there,  also  a  paper  by  Figueredo  to  appear)  which  gives  a  lower 
estimate  for  the  seventh  eigenvalue  (the  first  six  vanish).   This 
estimate  is  not  elementary  and  it  is  far  from  obvious  how  to 
extend  it  to  corresponding  finite  difference  equations  directly. 

Another  question  that  naturally  arises  concerns  estimates  of 
the  values  of  the  solution  at  Individual  points.   Such  estimates 
were  derived  from  the  maximum  principle  by  Gerschgorln  for  the 
Dlrichlet  problem  and  by  Batschelet  for  the  third  boundary  value 
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problem   [see  Forsythe,  Wasow] .   It  does  not  seem  possible  to 
establish  such  estimates  for  the  second  boundary  problem  by  using 
just  the  maximum  principle,  even  less  for  systems  where  no  simple 
maximum  principle  is  available. 

It  would  be  possible  to  obtain  such  estimates  by  the  method 
used  for  differential  equations  provided  the  finite  difference 
scheme  is  translation  invariant.   This  is  the  case  in  the  interior 
of  the  domain  provided  the  function  t^  =  (1/2)t]   +  (1/2)11   is 
employed.   Such  estimates  would  depend  on  derivatives  of  the  data; 
also  they  would  not  be  very  accurate.   Whether  such  estimates  could 
be  given  near  the  boundary  is  very  doubtful. 
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Appendix  II 

Remarks  about  the  Dirichlet  Problem 

In  this  appendix  we  shall  make  some  remarks  concerning  the  question, 
how  to  set  up  a  finite  difference  scheme  for  the  Dirichlet  problem  and  for 
problems  of  systems  in  which  all  boundary  values,    or  some  linear  combina- 
tions of  them,    are  prescribed. 

Suppose  the  boundary  happens  to  be  a  net  polygon  composed  of  "net 
segments",    i.  e.  ,    sides  or  diagonals  of  the  basic  rectangles  of  the  net;    and 
suppose  the  values  prescribed  on  the  boundary  are  linear  on  each  net  seg- 
ment.    Then  one  need  only  assign  the  prescribed  values  to  the  points  on  this 
polygon  and  otherwise  proceed  as  before,    provided  we  work  with  the  inter- 

polation  function    n      (or    n    ). 

~  I  II 

In  case  we  work  with  the  function    n  =  (l/2)r]    +  (l/2)r7       we  must  admit 

certain  external  points.     Let    P    P   P^P      o    be  the  end  point  of  a  net  rectangle 

and  suppose  that  the  diagonal    P    P^    is  part  of  the  boundary    B;      then  both 

points    P   ,     P  (one  of  which  being  exterior)  should  be  admitted.     In  addi- 

o         a+h 

tion  to  prescribing    u      and    u^    at     P      and    P^,     the  boundary  condition 

a  p  a  p 

(-^>  ^"'K    '    ^'^'K.^    -    "(l/2Ml/2)^ 

must  then  be  imposed.     Here    u      ,  ,     o    is  the  given  boundary 

value  at  the  midpoint  of  the  rectangle.     The  function  resulting 

from  interpolation  then  assumes  the  prescribed  piecewise  linear 

boundary  values.     A  simple  consideration  will  show  that  the  "positive 
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character  of  the  resulting  finite  difference  equation  will  not  be 
affected  by  the  midpoint  condition  (+) . 

For  the  case  that  the  boundary  cannot,  with  sufficient 
accuracy,  be  taken  as  a  net  polygon  the  following  procedure  is 
offered  for  consideration. 

A  net  polygon  B     is  chosen  which  approximates  the  actual 
boundary  (B  in  such  a  way  that  (£,   and  (B   can  be  mapped  into  each 
other  without  much  distortion.   Moreover  3   should  be  so  chosen 
that  at  most  two  vertices  of  a  net  triangle  belong  to  (B  . 

Let  P  be  the  point  on  CB  which  corresponds  to  the  netpolnt  P 
on  B  .   Then,  (B  will  consist  of  segments  joining  pairs  of  such 
points  P  and  these  will  correspond  to  the  sides  and  diagonals 
forming  the  polygon  (B  .   We  now  replace  the  segments  on  (B  by 
straight  segments  joining  the  same  points  P,  thus  forming  a 
polygon  (B-   We  then  prescribe  our  boundary  values  on  CB  instead  of 
©  by  interpolating  the  values  given  at  the  points  P  linearly.   The 
Dirichlet  problem  will  then  be  considered  not  in  the  region  Jl  but 
in  the  region  ^  enclosed  by  B. 

•X- 

Every  net  segment  connecting  a  net  point  P  not  on  B  with 
a  net  point  P   on  ©   is  now  replaced  by  a  straight  segment  connect- 

-x- 

ing  P  with  the  point  P  on  £  which  corresponds  to  P  .   Also  the 

■X-         * 

segment  connecting  two  points  P   on  ©  are  replaced  by  the  segment 
connecting  the  corresponding  points  on  (B-   Thus  a  modified  tri- 
angulation  of  the  plane  is  obtained. 
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The  external  points  which  must  be  admitted  in  case  we  work  with 
the  function  r;  should  not  be  changed  in  this  mapping. 

We  now  consider  functions  u  which  are  piecewise  linear  with 
respect  to  this  new  triangulation  and  which  assume  the  prescribed 
values  on  the  boundary  points  P.  Since  the  (modified)  prescribed 
boundary  values  on  the  (modified)  boundary  CB  are  piecewise  linear 
it  is  clear  that  any  function  as  described  assumes  the  prescribed 
boundary  values  on  <B  and  hence  is  admitted  to  the  variational 
problem  of  minimizln^^  the  integral 

D^  u  -  \\_   fudxdy 

where  (R  is  the  region  enclosed  by  3.   Variation  with  respect  to 
the  values  of  u  at  the  net  points  P  not  on  B  then  leads  to  a  finite 
difference  equation. 

In  order  to  interpret  the   nature  of  this  equation  we  find  it 
convenient  to  transform  the  region  <R  into  the  polygonal  region  3 
by  mapping  the  modified  triangles  into  the  corresponding  triangles 
of  the  original  net.   We  simply  map  the  points  P  into  the  corre- 
sponding points  P;  next  map  the  segments  connecting  points  P  with 
points  P  or  other  points  P"  linearly  into  the  corresponding  net 
segments;  and  finally  map  the  triangles  with  points  P  as  vertices 
linearly  into  the  corresponding  net  triangles. 

Since  the  mapping  of  each  of  these  triangles  is  linear,  the 

piecewise  linear  function  u  goes  over  into  a  piecewise  linear 

function,  but  the  derivatives  ^  u,  S  u  go  over  into  linear 

X    y 
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combinations  of  the  derivatives  of  u  as  functions  of  the  new 
variables.   The  coefficients  of  these  linear  combinations  are 
constant  in  each  triangle. 

In  any  case,  after  this  transformation  the  problem  has  been 
reduced  to  one  with  nonconstant  coefficients  in  a  polygonal  domain, 
thus  to  a  problem  of  the  type  treated  in  Section  4,  except 
that  now  the  coefficients  are  plecewlse  constant  and  not  contin- 
uous.  Formulas  (^  of  page  25  and  (^#)  with  ^J^    on  page   :.. ,  , 
27'  are  therefore  applicable  but  not  the  subsequent  comparison  of 
difference  and  differential  equations.   In  fact,  the  finite 
difference  expression  as  here  proposed  for  near  boundary  points 
does  not  go  over  into  the  differential  equation  if  the  mesh  width 
shrinks  to  zero. 

An  Important  question  is  whether  or  not  the  finite  difference 
equation  set  up  in  the  manner  described  is  of  positive  type.   This 
will  be  the  case  if  the  distortion  caused  by  our  transformation  is 
not  too  large.   The  question  whether  or  not  this  can  be  attained 
by  making  the  mesh  width  small  enough  (assuming  the  curvature  of 
the  boundary  to  be  finite)  must  for  the  present  be  left  open. 
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